##Generate data and maps for USDOL Prisoner Reentry Grant (April 2017)

##Load packages
library(tigris)
## 
## Attaching package: 'tigris'
## The following object is masked from 'package:graphics':
## 
##     plot
library(stringr)
library(xlsx)
## Loading required package: rJava
## Loading required package: xlsxjars
library(readxl)
library(leaflet)
library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(ggplot2)
library(ggmap)
library(RColorBrewer)
library(rgdal)
## rgdal: version: 1.2-5, (SVN revision 648)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/sschmid/Documents/R/R-3.3.3/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/sschmid/Documents/R/R-3.3.3/library/rgdal/proj
##  Linking to sp version: 1.2-4
##Set working directory to Shapefile folder
setwd("~/R Projects/USDOL_ReentryProject/Shapefiles")

##Read in saved shapefiles
####Michigan cities
cities <- readOGR(dsn = "Michigan_Cities", 
                  layer = "tl_2010_26_place00")
## OGR data source with driver: ESRI Shapefile 
## Source: "Michigan_Cities", layer: "tl_2010_26_place00"
## with 630 features
## It has 17 fields
## Integer64 fields read as strings:  ALAND00 AWATER00
####DPD precincts
precincts <- readOGR(dsn = "DPD Precincts", 
                     layer = "geo_export_7b552edd-0682-452c-b491-7e6a16b1c277")
## OGR data source with driver: ESRI Shapefile 
## Source: "DPD Precincts", layer: "geo_export_7b552edd-0682-452c-b491-7e6a16b1c277"
## with 12 features
## It has 12 fields
####Detroit neighborhoods
neighborhoods <- readOGR(dsn = "Detroit Neighborhoods", 
                         layer = "geo_export_9b67ef6c-9a29-4277-87eb-e2d8eafc5186")
## OGR data source with driver: ESRI Shapefile 
## Source: "Detroit Neighborhoods", layer: "geo_export_9b67ef6c-9a29-4277-87eb-e2d8eafc5186"
## with 201 features
## It has 1 fields
##Read in shapefiles from tigris package
####Southeast MI census tracts
tracts <- tracts(state="MI", county=c(99, 125, 163), cb=TRUE)
####Michigan counties
counties <- counties(state = "MI", cb = TRUE)
####ZCTA starting with "48"
zipcodes <- zctas(cb=TRUE, starts_with=c("48"))

##Set working directory to Raw Data folder
setwd("~/R Projects/USDOL_ReentryProject/Raw Data")

##Read in raw data files
####ACS 2014 5-year poverty data by census tract
tractfile <- read_excel("Poverty_CensusTracts_492017.xlsx")
####UWSEM school partner names and geocoordinates
schoolfile <- read.csv("UWSEM School Coordinates.csv", header=TRUE)
####DPD 2015 crime incidents
crimefile <- read.csv("Crime_2015.csv", header=TRUE)

##Set working directory to main project folder
setwd("~/R Projects/USDOL_ReentryProject/")


##Clean data and shapefiles
####Keep only SE Michigan counties in counties shapefile
counties <- counties[(counties$NAME == "Macomb" | 
                        counties$NAME == "Oakland" | 
                        counties$NAME == "Wayne"),]
####Remove non-violent crime data
crimefile <- crimefile[which(crimefile$CATEGORY == "AGGRAVATED ASSAULT" |
                               crimefile$CATEGORY == "HOMICIDE" |
                               crimefile$CATEGORY =="ROBBERY"),]
####Change neighborhood names in shapefile to match crime data
neighborhoods$name <- toupper(neighborhoods$name)
neighborhoods$name[neighborhoods$name=="BARTON-MCFARLAND"] <- 
  "BARTON MCFARLANE"
neighborhoods$name[neighborhoods$name=="SOUTHWEST DETROIT"] <- 
  "SOUTHWEST DETROIT / MEXICANTOWN"

##Merge data to shapefiles
####Merge ACS poverty rate data to census tract shapefile
merged <- geo_join(tracts, tractfile, "AFFGEOID", "AFFGEOID")
##Duplicate census tract shapefile and remove high poverty census tracts
mergedpoverty <- merged
mergedpoverty <- mergedpoverty[!(is.na(mergedpoverty$poverty)),]
mergedpoverty <- mergedpoverty[!(mergedpoverty$poverty > 30),]
####Merge crime data to census tract shapefile 
crimetract <- table(crimefile$CENSUSTRACT)
crimetract <- as.data.frame(crimetract)
crimetract$tract <-crimetract$Var1
mergedtract <- geo_join(merged, crimetract, "NAME", "tract")

##Calculate violent crime rate/50000 residents
mergedtract$crimerate <- (mergedtract$Freq/mergedtract$population)*50000

##Duplicate census tract crime file and remove high crime census tracts
mergedcrime <- mergedtract
mergedcrime <- mergedcrime[!(is.na(mergedcrime$crimerate)),]
mergedcrime <- mergedcrime[!(mergedcrime$crimerate > 840),]

##Merge crime data to Detroit neighborhood shapefile
crimeneighborhood <- table(crimefile$NEIGHBORHOOD)
crimeneighborhood <- as.data.frame(crimeneighborhood)
crimeneighborhood$name <-crimeneighborhood$Var1
mergedneighborhood <- geo_join(neighborhoods, crimeneighborhood, "name", "name")

##Merge crime data to DPD precincts shapefile
crimeprecinct <- table(crimefile$PRECINCT)
crimeprecinct <- as.data.frame(crimeprecinct)
crimeprecinct$PRECINCT <-crimeprecinct$Var1
####Clean precinct names in crime data file to match shapefile
crimeprecinct$PRECINCT <- as.character(crimeprecinct$PRECINCT)
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="2"] <- "02"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="3"] <- "03"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="4"] <- "04"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="5"] <- "05"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="6"] <- "06"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="7"] <- "07"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="8"] <- "08"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="9"] <- "09"
####Merge files
mergedprecinct <- geo_join(precincts, crimeprecinct, "name", "PRECINCT")

##Add population data to precinct shapefile - source: 2010 Census
####Create data frame with population data
population <- c("54499", "43701", "66398", "49435", "63701", "46342", 
                "86006", "74348", "56894", "56596", "75052")
precinct <- c("02", "03", "04", "05", "06", "07", "08", "09", "10", "11",
              "12")
precinctpop <- cbind(population, precinct)
precinctpop <- as.data.frame(precinctpop)
####Merge files
mergedprecinctcrime <- geo_join(mergedprecinct, precinctpop, "name", "precinct")
####Calculate violent crime rate/50000 residents for each precinct
mergedprecinctcrime$Freq <- as.numeric(mergedprecinctcrime$Freq)
mergedprecinctcrime$population <- 
  as.numeric(as.character(mergedprecinctcrime$population))
mergedprecinctcrime$crimerate <- 
  (mergedprecinctcrime$Freq/mergedprecinctcrime$population)*50000

##Create labels and popups
popupschool <-paste0(schoolfile$school)
zipcodelabel <- paste0(zipcodes$ZCTA5CE10)
popup1 <- paste0("Census Tract: ", merged$NAME, "<br>", 
                 "Population: ", merged$population, "<br>",
                 "Poverty Rate: ", round(merged$poverty, 2), "%")
popupcity <-paste0(cities$NAME00)
precinctlabel <- paste0("DPD Precinct ", precincts$name)
popuplabel <- paste0(neighborhoods$name)
popup2 <- paste0("Census Tract: ", mergedtract$NAME, "<br>",
                 "Population: ", mergedtract$population, "<br>",
                 "Violent Crimes: ", mergedtract$Freq, "<br>",
                 "Violent Crime Rate: ", round(mergedtract$crimerate,0))
popuphood <- paste0("Neighborhood: ", mergedneighborhood$name, "<br>",
                    "Violent Crimes: ", mergedneighborhood$Freq)
popupprecinct <- paste0("Precinct ", mergedprecinctcrime$name, "<br>",
                        "Violent Crime Rate: ", 
                        round(mergedprecinctcrime$crimerate, 0))

##Create color palettes for crime and poverty rates
pal <- colorNumeric(palette="Reds",
                    domain = merged$poverty, na.color="Gray")
qpal <- colorBin(palette="Reds",
                 domain = mergedtract$crimerate, na.color="Gray", 
                 bins=8)
palhood <- colorBin("Reds", domain=mergedneighborhood$Freq, bins=4)
palprecinct <- colorBin("Reds", domain=mergedprecinctcrime$crimerate)

##Create leaflet map
map1<-leaflet() %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  setView(lng=-83.170268, lat=42.373568, zoom=10) %>%
  addPolygons(data=counties, group="Census Tract Poverty Rates", 
              fillColor="Gray", color="black", 
              fillOpacity = 0, weight=5) %>%
  addPolygons(data = merged, group="Census Tract Poverty Rates", 
              fillColor = ~pal(poverty), 
              weight=0.3, fillOpacity = 0.7, color="Gray",
              popup = popup1) %>%
  addPolygons(data=counties, group="Census Tract Violent Crime Rates", 
              fillColor="Gray", color="black", 
              fillOpacity = 0, weight=5) %>%
  addPolygons(data = mergedtract, group="Census Tract Violent Crime Rates", 
              fillColor = ~qpal(crimerate), 
              weight=0.3, fillOpacity = 0.7, color="Gray",
              popup = popup2) %>%
  addPolygons(data=counties, group="DPD Precinct Violent Crime Rates", 
              fillColor="Gray", color="black", 
              fillOpacity = 0, weight=5) %>%
  addPolygons(data = mergedprecinctcrime, 
              group="DPD Precinct Violent Crime Rates", 
              fillColor = ~palprecinct(crimerate), 
              weight=0.3, fillOpacity = 0.7, color="Gray",
              popup = popupprecinct) %>%
  addLayersControl(baseGroups = c("Census Tract Poverty Rates", 
                                  "Census Tract Violent Crime Rates",
                                  "DPD Precinct Violent Crime Rates"),
                   overlayGroups = c("City Borders", 
                                     "Zip Codes",
                                     "DPD Precincts",
                                     "Detroit Neighborhoods",
                                     "High Poverty Census Tracts Only",
                                     "High Crime Census Tracts Only",
                                     "Returning Residents Center",
                                     "UWSEM Schools"),
                   position="bottomleft",
                   options = layersControlOptions(collapsed = FALSE)
  ) %>%
  addPolygons(data=cities, group="City Borders", fillColor="Gray",
              color="Black", fillOpacity=0, weight=2, dashArray="3",
              label=popupcity) %>%
  addPolygons(data=zipcodes, group="Zip Codes", fillColor="Gray",
              color="Black", fillOpacity=0, weight=2, dashArray="3",
              label=zipcodelabel) %>%
  addPolygons(data=precincts, group="DPD Precincts", fillColor="Gray",
              color="Black", fillOpacity=0, weight=2, dashArray="3",
              label=precinctlabel) %>%
  addPolygons(data=neighborhoods, group="Detroit Neighborhoods",
              fillColor="Gray", color="Black", fillOpacity=0, weight=2,
              dashArray="3",
              label=popuplabel) %>%
  addPolygons(data=mergedpoverty, group="High Poverty Census Tracts Only",
              fillColor="Gray", color="Gray", fillOpacity=0.7, weight=1) %>%
  addPolygons(data=mergedcrime, group="High Crime Census Tracts Only",
              fillColor="Gray", color="Gray", fillOpacity=0.7, weight=1) %>%
  addMarkers(lng=-83.170268, lat=42.373568, group="Returning Residents Center",
             label="Returning Residents Resource Center") %>%
  addMarkers(group="UWSEM Schools", lng=schoolfile$lng, lat=schoolfile$lat, 
             label=popupschool) %>%
  hideGroup("Zip Codes") %>%
  hideGroup("City Borders") %>%
  hideGroup("DPD Precincts") %>%
  hideGroup("Detroit Neighborhoods") %>%
  hideGroup("High Poverty Census Tracts Only") %>%
  hideGroup("High Crime Census Tracts Only") %>%
  hideGroup("UWSEM Schools")
## Warning in qpal(crimerate): Some values were outside the color scale and
## will be treated as NA
##View map
map1